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A key property of many-body localization, the localization of quantum particles in systems with 
both quenched disorder and interactions, is the area law entanglement of even highly excited eigen¬ 
states of many-body localized Hamiltonians. Matrix Product States (MPS) can be used to efficiently 
represent low entanglement (area law) wave functions in one dimension. An important application 
of MPS is the widely used Density Matrix Renormalization Group (DMRG) algorithm for finding 
ground states of one dimensional Hamiltonians. Here, we describe two algorithms, the Shift and 
Invert MPS (SIMPS) and excited state DMRG which finds highly-excited eigenstates of many-body 
localized Hamiltonians. Excited state DMRG uses a modified sweeping procedure to identify eigen¬ 
states whereas SIMPS is a shift-inverse procedure that applies the inverse of the shifted Hamiltonian 
to a MPS multiple times to project out the targeted eigenstate. To demonstrate the power of these 
methods we verify the breakdown of the Eigenstate Thermalization Hypothesis (ETH) in the many- 
body localized phase of the random field Heisenberg model, show the saturation of entanglement in 
the MBL phase and generate local excitations. 


Many-body localization (MBL) is a dynamical phase 
transition that happens at finite energy density for an 
isolated quantum system with interactions and quenched 
disorder CHZ], see also review [ 5 ]. Many-body localized 
phases are believed to have a number of properties includ¬ 
ing: (a) zero conductivity at finite temperature, (b) Pois¬ 
son statistics of their many-body eigenvalues, (c) eigen¬ 
states that obey the area law, (d) an extensive number 
of local integrals of motion in cases where there isn’t a 
mobility edge, and (e) eigenstates that fail to obey the 
eigenstate thermalization hypothesis (ETH) and hence 
MBL systems which fail to thermalize [HHH]- While the 
ground state wave-function is the key quantity for identi¬ 
fying equilibrium quantum phases, the finite-energy den¬ 
sity eigenstates are the analogous identifying feature for 
dynamical phases. Analytical arguments for the exis¬ 
tence of MBL phases and phase transitions rely primarily 
on analyzing the highly excited eigenstates using either 
diagrammatic re-summation [iHain] or the real space 
renormalization group From the numerical per¬ 

spective, getting interior eigenvalues (and corresponding 
eigenstates) using exact diagonalization requires both an 
exponentially large amount of computer time and mem¬ 
ory, and as a result these are limited to rather small (at 
most 22-site m) systems m uMi]. Time evolution 
DMRG studies, on the other hand, are limited by the 
logarithmic-in-time growth of entanglement entropy that 
occurs in the localized phase . It has been recently 

shown that the entire spectrum of eigenstates of a MBL 
system can be efficiently described by a matrix product 
operator (MPO) [551 [55]; the best available algorithm for 
optimizing this MPO m captures features of the entire 
spectrum but doesn’t describe individual states to high 
fidelity. While the bulk of MBL research has been the¬ 
oretical, there has been considerable recent progress on 
the experimental side [551 [ 55 ]. For the first time, clear 


signs of the MBL transition were observed in an ultracold 
atom experiment [30]. 

In this work, we take advantage of the fact that eigen¬ 
states of many-body localized systems obey the area law 
and can therefore be efficiently represented as a Matrix 
Product State (MPS) [51I3I1I31] to develop a set of numer¬ 
ical algorithm for generating a MPS representations of 
these eigenstates. We use the eigenstates constructed us¬ 
ing these algorithms to test the basic properties of MBL, 
the break-down of ETH, the saturation of the entangle¬ 
ment entropy, and the existence of a large number of 
local excitations, in the regime of large one-dimensional 
systems that were previously inaccessible due to the lim¬ 
itations of exact diagonalization. 

The algorithms we develop fall into two broad classes: 
In class (1), we modify the DMRG sweeping proce¬ 
dure [33] to pick one of the excited eigenstate of the ef¬ 
fective Hamiltonian at each step and hence arrive at an 
excited eigenstate of the full Hamiltonian. We call this al¬ 
gorithm Excited State-DMRG (ES-DMRG). In class (2), 
we target a specific energy A by repeatedly applying the 
operator {H — A)“^ to the state vector. We call this 
algorithm shift-and-invert MPS (SIMPS). SIMPS is en¬ 
abled by the inverse-DMRG algorithm, that we develop 
here, for efficiently computing the single application of 
the propagator {H — A)“^ to a MPS. 

Throughout this manuscript, we use the following one¬ 
dimensional Hamiltonian 

H = Y,Si-S,+i+Y.h,S^ ( 1 ) 

i i 

where hi are sampled uniformly from [—IT, IT], which is 
known to have the entire many-body spectrum become 
localized for IT > 3.5 mm- After developing our al¬ 
gorithms and showing that they generate eigenstates, we 
use the eigenstates to test a number of physical proper- 


2 


ties. First, we verify that in the many-body localized 
phase the eigenstate thermalization hypothesis (ETH) 
breaks down. We explicitly demonstrate this using eigen¬ 
states in a tight energy window that we find using SIMPS 
for L = 30 chains. The second key feature that we ver¬ 
ify is the saturation of the entanglement entropy with 
system size. We verify that this saturation is occurring 
using eigenstates obtained via SIMPS and ES-DMRG for 
L = 30,40. The third key feature of MBL phases, that we 
investigate are the local excitations. We show how, given 
an eigenstate, we can construct approximate eigenstates 
that differ from a reference eigenstate locally. While this 
procedure does not have enough power to construct a 
full set of mutually orthogonal local integrals of motion, 
it does produce a large number of of “pretty good” eigen¬ 
states that can potentially be used to investigate the local 
integrals of motion. 

Overcoming the exponential gap: As a physical 
system approaches the thermodynamic limit, the inter¬ 
level spacing decreases exponentially. Therefore expo¬ 
nentially many states exist within any fixed energy win¬ 
dow {E — 5,E + 5). The fundamental difficulty in cap¬ 
turing eigenstates, then, is resolving a single eigenstate 
from a large superposition of eigenstates taken from this 
window. 



FIG. 1. Average mid-bond entanglement entropy of the wave- 
function = cos(a)'i/)i -I- sin(a)'!^i+i as a function of the 
mixing angle a, where tpi and 'ipi+i are consecutive eigenstates 
taken from the middle of the spectrum of the Hamiltonian 
in Eq. Q. Notice that pure eigenstates are local minima 
of the mean entanglement entropy even when approaching 
no disorder. The average is over 1% of eigenstates on a 14 
site chain; a single realization of disorder was used for each 
disorder strength. 

While it is well known that MBL eigenstates (MBLE) 
obey area laws [9] , we show that MBLE are local minima 
in entanglement with respect to linear superpositions of 
eigenstates nearby in energy (see Eig. [^. Conceptually 
this happens because nearby MBLE are typically local¬ 
ized in different parts of the Hilbert space and building 
a superposition of them produces a cat state with larger 
entanglement. This feature gives an additional metric 
from which we can identify the MBLE states from the 
nearly degenerate subspace: we want to select the min¬ 


imally entangled states or equivalently MBS with low 
bond-dimension in this space. This realization allows us 
to develop successful numerical algorithms for studying 
these eigenstates. 

Excited state DMRG: The first of these algorithms, 
while conceptually interesting, is simple to describe and 
has the advantage of requiring minimal modifications to 
a current DMRG code. A standard DMRG code sweeps 
over sites; when working on site z, an effective Hamilto¬ 
nian Hi is generated which involves tracing over the other 
auxiliary and physical degrees of freedom. This effective 
Hamiltonian Hi is then diagonalized and the parameters 
on site i are replaced with the ground state of Hi. In our 
new approach, instead of considering only the ground 
state of Hi we consider all its eigenvalues (of which there 
are replacing the parameters of site i with the one of 
these eigenvalues. In particular, we select the eigenvalue 
with energy closest to the energy of the current MBS to 
minimize the amount of change in the state as our algo¬ 
rithm progresses. We then follow the typical approach of 
sweeping back and forth through all the sites. As in nor¬ 
mal DMRG, we find that a proper starting configuration 
and bond-dimension protocol can enhance the efficiency 
of the algorithm keeping it from being stuck and allow¬ 
ing it to more widely sample excited states. In partic¬ 
ular, we typically start with the algorithm in a product 
state in the Sz basis of bond-dimension two and slowly 
ramp up the bond-dimension increasing it by one every 
few sweeps. This ensures we find low bond-dimension 
states at the energy at which the algorithm converges. 
Erom these sweeps, we take the state we find with low¬ 
est variance. This algorithm scales as typical DMRG. 
While currently we are exactly diagonalizing the effec¬ 
tive Hamiltonian, a standard shift-and-invert procedure 
would allow the eigenstate of the effective Hamiltonian 
to be found without a full diagonalization. While this 
approach is powerful, it has some undesirable properties. 
In particular, there is no clean way to target a particular 
energy. 

Shift and Invert MPS (SIMPS): Our goal is to 
target a particular energy A. The simplest approach is 
to use standard DMRG techniques to find the ground 
state of {H — Xf' instead of the ground state of H. Un¬ 
fortunately this decreases an already exponentially small 
gap in nearby eigenvalues by squaring it. As techniques 
such as ITEBD require propagating to imaginary time 
larger then this inverse gap, they are unfeasible; we have 
found sweeping methods are able to converge to some 
eigenstate but appear to find one further from the en¬ 
ergy target than the SIMPS approach described below. 

Exact diagonalization is also plagued by the small gap 
coming from working with (H—E)'^. Instead of using this 
propagator the largest MBLE which have been generated 
via ED are N = 22 [15] and use the shift-and-invert tech¬ 
nique; this technique repeatedly applies {H — A)“^ to 
a random state converging to an eigenstate with eigen- 
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value close to A. This technique overcomes the problem 
of nearly degenerate eigenstates by inverting the spec¬ 
trum. We develop an analogous approach in the MPS 
language called SIMPS. Starting with a random MPS, 
we iteratively apply {H — A)“^ until we reach an MPS 
that well approximates an eigenstate close to A. The con¬ 
vergence speed of this method is geometric in the limit of 
large bond dimension, and mostly controlled by the ratio 
of the second dominant eigenvalue to the first dominant 
eigenvalue present in the state 


bd 

to 

1 
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1 
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E 2 — X 


( 2 ) 


where Ei is the eigenvalue of H closest to A and E 2 is 
the second closest to A. As an example, given an MPS 
IV') = a\Ei) + h\E 2 ) + ■ • • (where • • • means other states 
far away from A), the MPS \ip) = {H — X)~^\tjj) has an 
energy of 


(E) 


1 , M 

lap- 


(3) 


Therefore, the energy of the MPS decays exponentially in 
N. To increase the convergence speed, one can fine tune 
A such that p is much smaller than 1, by fitting the energy 
trace of (E) vs. N accumulated during a particular run 
to Eq. See Supplementary Figures for a prototypical 
example of SIMPS converging and an example of this 
decay with a poorly chosen value of A. 

Inverse DMRG: The key to SIMPS is the develop¬ 
ment of an inverse-DMRG approach which variationally 
applies the inverse of an MPO to an MPS. Just like in¬ 
verting a non-degenerate matrix A directly and applying 
it to a vector b is more difficult than solving Ax = b, 
inverting an MPO directly suffers two major problems - 
(1) it lacks an efficient algorithm; (2) it may potentially 
require very large bond dimension for the inverted MPO. 
The more efficient alternative is to construct a MPS \ip) 
which variationally approaches 0“^|V’), for a given MPO 
O and MPS |V’)- To this end, we propose to build the 
MPS \ip) by minimizing IIOI-/?) — IV’)IP; 

^||O|<p)-|V’)|P=0, (4) 


where i stands for a particular site and cr is a label for 
degrees of freedom on site i. The above equation leads 
to the following variational prescription 

(5) 
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FIG. 2. Pictorial representation for optimizing ||0|(/3) — IV-)!p. 
Here as an example, we consider optimizing the third site 
(colored orange) for a L = 4 chain. The MPO O (pink) and 
the MPS IV) (cyan) are given. To update the orange block, 
we solve a linear equation, treating the orange block as an 
unknown vector x\ the other part of the network on the left 
hand side, after being contracted, amounts to a symmetric 
semi positive-definite matrix A, and the matrices on the right 
hand side of the equation becomes a known vector b. The 
entire MPS is optimized by sweeping. Given an cutoff e, the 
sweeping can be stopped by checking if |1 — ((p|0^|V’)| < £, 
which can be done without additional computational effort 
since can be calculated by contracting the solved 

orange block with the rest of the network on the right. 


MPS \ip), one needs to sweep the chain back and forth as 
in standard DMRG. 

Precision of this method is controlled by the bond di¬ 
mension of the MPS |(^), and the number of sweeps. The 
computationally dominate piece is in solving the dense 
system of linear equations. Substantial speedup can be 
gained by using conserved quantum numbers (QN) (al¬ 
though we currently don’t use this) which reduce the size 
of equations. The computational cost of inverse-DMRG 
scales as 0{L{pM‘^)^) for using a direct solver, and op¬ 
timistically when using an iterative solver, 

where p is the number of physical degrees of freedom per 
site, and M is the maximum bond dimension of the MPS. 
We currently do direct solves up to M = 60. 

Generating Eigenstates: We start by producing 
eigenstates at length L = 10 such that they can be com¬ 
pared with ED. Here, for both algorithms, we consider a 
fixed disorder configuration with W = 8 and artificially 
limit bond dimension to M = 12. Since ES-DMRG can’t 
target particular energies, we run the algorithm many 
times verifying that the energies it finds match those of 
the true Hamiltonian (see Fig. [^. In SIMPS, we have a 
tunable shift parameter A and focus on the energy win¬ 
dow of A G [—0.1,0.!], within which there are 12 eigen¬ 
states as shown by ED. When running at a limited bond 
dimension with fixed initial conditions, the SIMPS algo¬ 
rithm does not always hit the targeted eigenstate. Re¬ 
markably, it does find an eigenstate near the desired en¬ 
ergy with high fidelity. In fact, if we consider 

P(A) =min[l-|(MPS'(A)|£;D(i))|] (6) 

i 


which is described pictorially in Fig. To update the 
matrices at a particular site, one needs to solve a lin¬ 
ear equation problem that involves a dense, symmetric, 
and semi positive-definite matrix. To optimize the entire 


then for L = 10, M = 12 we never see a situation where 
the P{X) > 10-^ Fig. I shows the eigenvalues identi¬ 
fied during a A sweep for different initial conditions and 
bond-dimensions of SIMPS. Notice that changing the ini- 
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FIG. 3. Convergence to eigenstates 503-514. Magenta lines 
indicate the twelve eigenvalues in this range. SIMPS: Blue 
and yellow bars indicate the basins of convergence, as a func¬ 
tion of parameter A, found for two random MPS initial con¬ 
ditions. Yellow and blue dots indicate whether sweeps of the 
two respective initial conditions find a given eigenvalue. ES- 
DMRG: Green dots indicate whether ES-DMRG found the 
eigenvalue after being run with 6144 random product state 
initial conditions. 


tial conditions allows targeting of eigenstates that were 
inaccessible with the original initial condition. Hence, 
using just two initial conditions it is possible to target 10 
out of the 12 eigenstates. 

Another useful metric to quantify the quality of eigen¬ 
states, that can be applied to system sizes greater then 
those accessible to ED, is the standard deviation of the 
energy a = a/ {H)'^ — Note that while a true eigen¬ 

state has cr = 0, numerically calculated standard devia¬ 
tions are limited by machine precision. Unsurprisingly, 
the high fidelity we see on short chains {L = 10, M = 12) 
corresponds to (logj^Q a) « —6.1 which is beyond machine 
precision. 

As both SIMPS and ES-DMRG algorithms stochasti¬ 
cally converge to eigenstates, one has to check whether 
the states that our algorithms converge to are typical or 
biased. Specifically, we want to verify there isn’t a bias 
toward eigenstates with lower entanglement. To verify 
this is not occurring we compare histograms of entan¬ 
glement entropies for eigenstates generated using ED to 
those generated using SIMPS and ES-DMRG at artifi¬ 
cially reduced bond-dimension M=12 (see Fig. [^. From 
the comparison we observe that both SIMPS and ES- 
DMRG are sampling states at a given entanglement with 
the same frequency as ED and hence there is no system¬ 
atic bias. 

We find that, using SIMPS at M = 60, we can obtain 
eigenstates of chains of length L — 30,40 and IT = 8 to 
average values of (cr) below machine precision and with a 
distribution that is largely peaked at machine precision 
(see fig. . These tests are conducted for the worst case 
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FIG. 4. Values of (logjQ((T)) for various algorithms and pa¬ 
rameters at different disorder strengths W. Blue line (L = 30) 
and orange line {L — 40) are SIMPS runs at M = 60 with A 
tuned exactly to the middle of the spectrum for each disor¬ 
der realization, where the lowest and highest eigenvalues are 
obtained via DMRG and the disorder conhgurations are the 
same scaled pattern at different W for L = 30 and L — 40, 
respectively. Values of IT < 5 are not necessarily reliably 
converged. Green point is ES-DMRG at M = 20, L — 30 and 
pink point is ES-DMRG at M = 20, L — 100. The dark green 
point is SIMPS at M = 20, L = 30 showing that ES-DMRG 
and SIMPS converge to similar accuracy at the same bond 
dimension and chain length. Insets are distributions of a at 
the respective parameters. 


scenario - we target A exactly to the middle of the en¬ 
ergy spectra (lowest and highest energies calculated by 
DMRG) where the many-body density of states is near 
maximum. Running ES-DMRG at M = 20 for chains 
with {L — 30, W = 8} and {L = 100, W = 10}, we find 
approximate eigenstates with (log]^o('^)) somewhat larger 
than SIMPS. 

Of course, for long chains, the inter-level spacing is 
going to be smaller then the a which can be resolved by 
machine precision and so it is useful to have an additional 
measure to check that we are not converging to a linear 
combination of eigenstates. It is expected that when tun¬ 
ing A with a small increment, two successive MPSs built 
by the method should either be the same eigenstate or 
be orthogonal to each other. This means the hdelity 

F{X,) = \{MPS{X,)\MPSiX,+i))\ (7) 

should alternate between zero and one. If instead, each 
independent simulation at A produced a different lin¬ 
ear superposition of nearby eigenstates we would expect 
F{Xi) to vary continuously between zero and one. The 
inset of Fig.j^shows eight distinct eigenstates with F{Xi) 
alternating between 1 (the same state) or 0 (an orthogo¬ 
nal state). 

Many-body localization in large systems: We 

now use the tools that we have developed to test, in a 
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FIG. 5. A comparison of mid-bond entanglement en¬ 
tropy histograms generated using ED, SIMPS, and ES-DMRG 
eigenstates. To generate the histograms, we fixed the disor¬ 
der strength at IT = 8, chain length at L = 10, and the 
maximum bond dimension used in SIMPS and ES-DMRG at 
M — 12. Next, we generated 101 disorder realizations for ED 
and SIMPS (ES-DMRG data is for 320 different disorder real¬ 
izations). Using these we construct all 103,424 eigenstates us¬ 
ing ED, 14,017 eigenstates using SIMPS, and 951 eigenstates 
using ES-DMRG. Error bars correspond to points rii ± y/nl, 
where n; is the number of entanglement entropies in the i- 
th bin. The error bars should be treated as a lower bound 
as they do not take into account the fact that for the same 
disorder configuration the mid-bond entanglement entropy of 
different eigenstates is correlated, and hence the number of 
independent samples is lower than what is expected by the 
naive estimate above. 





FIG. 6. Energy vs A for L = 30, IT = 10. 29 MPSs obtained 
using SIMPS with M = 40 are displayed on this plot, with 
8 of them having distinct energies. Inset: Fidelity function 
F{\) as define in Eq. 0. F{X) also shows 8 distinct states, 
and matches exactly with the energy curve. 

previously inaccessible regime of long chains, three key 
properties of MBL matter: failure to thermalize, low en¬ 
tanglement entropy of highly excited eigenstates, and a 
large number of local excitations. 

The eigenstate thermalization hypothesis (ETH) states 
that in thermalized quantum systems nearby eigenstates 
have very similar expectation values of local observ¬ 
ables We test this by computing the expecta¬ 

tion value of a representative local observable (S^) for 


FIG. 7. Left: Site value of (S^) as a function of energy 
at L = 30 for a fixed disorder configuration with disorder 
strength IT = 10. Eigenstates are obtained using SIMPS, 
with M = 40. Note the small energy range. The energy scan 
was performed with 200 A evenly spaced in [—0.1,0.!]. In 
the Sj = —3,—2,—1,0 sectors, a total of 74 distinct MPSs 
were found. Right: Site value of (5f) for L = 16, IT = 0.8 
generated using ED. 


two chains: one with L = 16 at IT = 0.8 (expected to 
obey ETH) and another with L = 30 at IT = 8 (expected 
to violate ETH). Specifically, we choose eigenstates from 
a small energy window in the range [-0.1,0.!]. We also 
filter the eigenstates by total S'J = S'f. For the case 
IT = 8, SIMPS discovered 74 eigenstates that were close 
in energy (but not consecutive). As shown in Fig. 
the corresponding (51) vary wildly confirming the vio¬ 
lation of ETH. On the other hand, for IT = 0.8 ETH 
prevails - (5|) lies in a narrow band, with the position 
of the band depending on the global conserved quantity 
5J [see Fig.g. 



FIG. 8. Dependence of mid-bond entanglement entropy S on 
the disorder strength for various chain lengths. 50 disorder 
realizations were used for L = 14 and 100 for L — 30 and L = 
40. Open markers indicate points at low disorder strengths 
for which SIMPS failed to obtain states with low standard 
deviation of the energy. For SIMPS we set M = 60 and tuned 
A to the exact middle of the eigenspectrum for each disorder 
realization. 
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The second key property that we investigate is the 
mid-bond entanglement entropy of highly excited eigen¬ 
states. In the ergodic phase these eigenstates have exten¬ 
sive mid-bond entanglement entropy while in the MBL 
phase the mid-bond entanglement entropy should satu¬ 
rate as a function of system size. In Fig. we plot the 
mean mid-bond entanglement entropy as a function of 
disorder strength for chains of length L = 14 (ED) and 
L = 30,40 (SIMPS). We observe that for IF >4 the mid¬ 
bond entanglement entropy is essentially independent of 
systems size, thus confirming the predicted saturation in 
the MBL phase. On the other hand for IF < 4 the ED 
and SIMPS data strongly disagree. We attribute this 
disagreement to the failure of SIMPS to hnd high qual¬ 
ity eigenstates in this near-ergodic and ergodic regime as 
demonstrated in Fig. 
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FIG. 9. Spectrum of 62 excitations (colored bars) produced 
by local modifications of a reference MPS eigenstate (long 
magenta line). The horizontal axis marks the site on which 
the MPS is being modified. Bars of the same color correspond 
to the same excitation as determined by their energies. \L = 
30, IF = 10, M = 15, excited eigenstates all have cr < 0.0005]. 

The third property that we investigate is the existence 
of a large number of local excitations of the system. In 
Fig. [^we show a spectrum of these local excitations con¬ 
structed with respect to an eigenstate at E = —15.002. 
To construct a local excitation on a given site, we take 
the following steps. (1) We bring the reference state to 
canonical form with respect to the target site. (2) We 
obtain the effective Hamiltonian on the target site. (3) 
We construct an MPS for all eigenstates of the effective 
Hamiltonian on the target site. (4) Most of the MPSs 
produced in step (3) are superpositions of a large num¬ 
ber of local excitations, and hence are not necessarily 
good eigenstates. Therefore, we retain only the MPSs 
that have a low a. Surprisingly, even after the filtering 
step, we find a large number of local excitations (there 
are 62 distinct ones in Fig.[^. 

Outlook: In our manuscript, we have presented two 
algorithms ES-DMRG and SIMPS for Ending matrix 


product state representations of inner eigenstates. We 
demonstrated the functionality and efficiency of both 
algorithms by showing that they can be used to pro¬ 
duce good variational MPS approximations of the ex¬ 
cited states of a fully many-body localized Hamiltonian 
in one dimension. We have also used these algorithms to 
demonstrate three key attributes of many-body localiza¬ 
tion, the breakdown of the eigenstate thermalization hy¬ 
pothesis, the saturation of the entanglement entropy, and 
the presence of a large number of local excitations, in a 
previously inaccessible regime of long chains, L = 30,40. 

We are particularly excited about two potential conse¬ 
quences of our work for many-body localization. First, 
as demonstrated in Fig. the algorithms that we have 
developed seem to work right up to the localization- 
delocalization transition. We speculate that with further 
improvements these algorithms can be used to probe the 
transition from the localized side m- Second, it seems 
possible that by collating and orthogonalizing the “local 
excitation” data of Fig. [^one can construct and charac¬ 
terize the local integrals of motion that are central the 
many-body localization phenomenology. 

We believe that the methods that we have developed 
are not limited to the study of strongly-localized mat¬ 
ter. Specifically, these methods have the potential to be 
significantly better for finding low-lying excited states 
of conventional Hamiltonians as compared to current 
state of the art methods - an extremely important prob¬ 
lem both in quantum chemistry and condensed matter 
physics. Moreover, we suspect that further improvements 
can be made to our algorithms by applying more sophisti¬ 
cated diagonalization methods to matrix product states. 

Note Added: During the preparation of this 
manuscript we became aware of two independent works 
that also developed an algorithm similar to ES-DMRG 
for finding excited eigenstates in MBL systems [35j and 
molecular systems [55] . 
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Number of power iterations 


FIG. 10. Convergence behaviors of one SIMPS run for L = 30, 
IF = 8 and M — 40. Initial MPS, at iteration 0, is a random, 
normalized MPS. Top to bottom: energy E, log of standard 
deviation logjQ((7), mid-bond von Neumann entanglement en¬ 
tropy S. In inverse-DMRG, when starting the sweeping for 
\ip) = {H — A)“'^|t/)), one could (a) start with \ip) the same as 
IV'); (b) start with random \ ip). The second approach can help 
to overcome local minima during inverse-DMRG sweeping. In 
the plot above, \ip) was set to a random MPS at iteration 1 
and 20. 
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